Order-Parameter Fluctuations in Ising Spin Glasses at Low Temperatures 
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We present a numerical study of the order-parameter fluctuations for Ising spin glasses in three and 
four dimensions at very low temperatures and without an external field. Accurate measurements of 
two previously introduced parameters, A and G, show that the order parameter is not self-averaging, 
consistent with a zero-temperature thermal exponent value 6' ~ 0, and confirm the validity of the 
relation G = 1/3 in the thermodynamic limit in the whole low-temperature phase, as predicted by 
stochastic stability arguments. 

PACS numbers: 75.40.Mg, 75.10.Nr 



I. INTRODUCTION 

Understanding the low-temperature physics of short- 
range spin glassesi remains a major unsolved problem. 
Much of the current debate concentrates on the equilib- 
rium thermodynamics of the Edwards-Anderson model 
with Ising spins (EAI model), the canonical short-range 
spin glass. Since analytical approaches pose formidable 
difficulties, the problem is often studied numerically. 
However, the existence of large barriers between low- 
energy configurations has limited so far numerical cal- 
culations to small systems sizes, from which it is hard to 
draw definite conclusions on the large- volume limit. 

Two main issues have been addressed in many numer- 
ical studies: the existence and character of the finite- 
temperature spin-glass transition, and the nature of the 
low-temperature spin-glass phase. A central quantity 
of interest in the description of the spin-glass phase is 
the scaling exponent 9' governing the typical energy of 
the lowest-lying excitations with linear size of order I, 
which is assumed to scale as E ^ . In general, 9' 
may be distinctSi^ from the stiffness exponent 9 mea- 
sured in domain-wall computations^i^i^, and the stabil- 
ity of the spin-glass phase requires to have 6*' > 0. In 
a "many-state" picture^, such as the replica-symmetry- 
breaking picture inspired by mean-field theorjiSiifi, one 
has 9' = 0, hence there are excitations whose energy 
remains finite (of the order of the coupling strength be- 
tween two spins) even as their length scale diverges. In a 
"two-state" picture, such as the droplet modelii*iS, one 
has 9' > 0, hence the energy of large-scale excitations 
diverges with their size. In this case the identity 9' = 9 
is often assumed. 

Both the spin-glass transition and the ordered phase 
have been usually investigated numerically by com- 
puting sample-averaged quantities such as the Binder 
cumulant^-i or the distribution of the order parameter 
(OP) and related observables. Recentl y^^'^^ , it was ob- 
served that useful information on both issues can be 
drawn from the sample-to-sample fluctuations of the OP. 



In particular, two dimensionless measures of the OP fluc- 
tuations were consideredi^*i^: A, the normalized fluc- 
tuation of the spin-glass susceptibility, and G, a ratio 
between two cumulants of the OP distribution. These 
two parameters are related to the Binder cumulant, B, 
via the relation B ~ 1 — A/{2G). For a model without 
time-reversal symmetry (TRS), A, G, and B are given by 
Eqs. 0, ®j and l@J below. The parameter G serves as 
a good indicator of the existence of phase transitions — 
in systems lacking TRS (for which B is generally a bad 
indicator), as recently shown for several systems, includ- 
ing Migdal-Kadanoff spin glassesi^, RNA folding mod- 
els^^, chiral spin systemsi^, and mean-field models such 
as the SK model (with and without a magnetic field), 
the inflnite-range p-spin model^^, and the inflnite-range 
Potts model2£. The parameter A has also been studied 
before for random diluted models at criticality^--^- . 

In this paper, we investigate the OP fluctuations in the 
EAI model with Gaussian couplings in three and four di- 
mensions by low-temperature Monte Carlo simulations. 
We study the case with no external field, which satis- 
fies TRS. In three dimensions (3D), numerical data are 
available in the literature for A in the high-temperature 
phase'^^ and for G near the critical point^'*, for the "± J" 
coupling distribution. In four dimensions (4D), G was 
measured at moderately low temperatures, also for the 
± J distributioni^. Here, we study much lower tempera- 
tures than in these studies, in order to reduce crossover 
effects associated to the critical point^^t^^, which compli- 
cate the interpretation of the numerical data at higher 
temperatures. 

A summary of our results is as follows. First, we esti- 
mate 9' from the system-size dependence of A, finding, 
for the system sizes we could reach, a small value of 9' 
incompatible with the accepted values of the domain- wall 
exponent [9 ~ 0.2 in 3D^ and 9 ~ 0.7 in 40^), and 
compatible with zero. This agrees with recent determina- 
tions of 9' from ground-state perturbation methodsSiSiSi 
and from low-temperature measurements of the OP dis- 
tributiorkifliSSiSSiSS (which all consider sample-averaged 
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quantities), and supports a picture of the spin-glass phase 
characterized by two distinct exponents, 9 > and 
61' = 0. The resuh 9' = implies that the OP is not 
self-averaging in the thermodynamic limit. 

Second, we find good evidence that the identity G = 
1/3 holds in the whole spin-glass phase in the thermody- 
namic limit, confirming the validity of sum rules proposed 
by Guerra^*' and first derived for the SK model, which 
follow from the property of "replica equivalence"i^ii2S. 

Third, we find that A and G allow to locate the spin- 
glass transition reasonably well although, as expected 
due to TRS and as previously numerically observed^, B 
provides a better determination (a much more accurate 
determination is provided by the correlation length24, 
which we do not investigate here). 

We do not study in this paper the surface fractal di- 
mension ds of the excitations, which is the other ex- 
ponent, besides 9', characterizing the spin-glass phase 
(in particular, = in the standard replica-symmetry- 
breaking picturei^, > in the droplet model, while the 
"TNT" picture^i^ predicts the "mixed" behavior dg > 0, 
9' ^0). 

The rest of the paper is organized as follows. In Sec. II 
we introduce the different models and observables stud- 
ied, and discuss the theoretical predictions for these ob- 
servables. In Sec. Ill we present and analyze our numer- 
ical results for the quantities A, G and B. Finally, in 
Sec. IV we summarize our conclusions. 



II. MODELS, OBSERVABLES, AND 
THEORETICAL PREDICTIONS 

We study the EAI model defined by the Hamiltonian 

5, = ±1 , (1) 

where Ising spins Si sit on a (hyper-)cubic lattice in 
D dimensions with linear size L and periodic boundary 
conditions in all directions. The couplings Jij are drawn 
from a Gaussian distribution of zero mean and unit vari- 
ance. We consider two different models: (i) the case with 
interactions restricted to nearest neighbors (referred 
to as NN model) in D = 3 and 4; (ii) the case with in- 
teractions restricted to nearest, next-nearest, and next- 
next-nearest neighbors (referred to as NNN model) in 
D = 3, which has a coordination number z = 26. 

The NN model has been extensively studied and is 
known to display a finite-temperature continuous spin- 
glass transition for D > 3. Recent estimates of the crit- 
ical temperature for Gaussian-distributed couplings are 
Tc = 0.95 ± 0.04 in 3DS and = 1.80 ± 0.03 in 4D2i. 

The NNN model has been much less studied. In 
Ref. [2J the 3D case with ± J couplings was considered, 
but no conclusive evidence of a finite-temperature transi- 
tion was obtained, the data being compatible with both 
Tc ~ 3.27 and a zero-temperature singularity. Inciden- 
tally, in the NNN case we do not expect a large difference 



in Tc between binary and Gaussian distributions, due 
to the large coordination number. For the NNN model, 
we did not consider temperatures as low as for the NN 
model, but focused on the phase transition region. An- 
other of the results of this paper is a convincing evidence 
that indeed a finite-temperature transition exists in this 
model in 3D. 

We measure A, G, and _B as a function of temperature 
T and size L using the following definitions: 

AiT,L) = Ml^JZl! (2) 

G(r,L) = (3) 

B{T,L) = 1(3-^1 (4) 

where q = SfS\ is the spin overlap of two inde- 

pendent systems Sf and S][ with the same random cou- 
plings, and (...) and (...) stand for thermal and disorder 
averages, respectively. The OP for a given realization 
of the disorder is (g^), therefore A is nothing but the 
normalized sample-to-sample variance of the OP. 

In the paramagnetic phase, T > Tc, and for sufficiently 
large L so that L » ^ (where ^ is the correlation length), 
the OP follows a Gaussian distribution and all three pa- 
rameters vanish as as follows from the central 
limit theorem. Following the terminology of Wiseman 
and Domany^^, this means that the OP is strongly self- 
averaging. 

At T — Tc, the correlation length diverges and the 
central limit theorem cannot be applied. For strongly 
disordered systems such as spin glasses it is known that 
the OP is not self-averaging at criticalit y^^' ' ^^ , namely A 
tends to a finite value in the thermodynamic limit. If 
A is finite then clearly G must be finite, and standard 
renormalization-group arguments show that B is also fi- 
nite at Tc- Since B and G are dimensionless and mono- 
tonic in T, in plots of these quantities as a function of 
T, the curves for different values of L must all cross at 
T — Tc, and one can use this to determine T^. From 
standard finite-size scaling, one can then determine the 
critical exponent v. Since much work has been devoted 
to measuring v from the sta ndard observables (see for 
instance Refs. I24'33'35' 3(l37^ . and we are primarily in- 
terested in the low-temperature phase here, we will not 
attempt a precise determination. 

In the spin-glass phase, T < Tc, A is expected^Si^S to 
vanish linearly with T according to the scaling law 

AiT,L)r.TL-'>' , (5) 

where 9' is the exponent discussed in the Introduction. 
This law holds under two hypothesis (both satisfied in 
the case of continuous couplings studied here): (i) the 
ground state is unique; (ii) the probability distribution 
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of the energy of the lowest-lying excitations has finite 
weight at zero energji22i^. 

From the above scaling law we see that if 6*' > 0, then 
A vanishes for L — > c«, namely the OP is weakly self- 
averaging (where "weakly" indicates^^ that the OP fluc- 
tuations vanish more slowly than a consequence of 
the inequality 9' < d). This situation is encountered in 
the droplet modeliiii^, as discussed in the Introduction, 
and also in mean-field models with a marginally stable 
replica-symmetric solution at low temperatures (such as 
the spherical SK model'"^). If 6*' = 0, as in a "many- 
state" picture, A remains finite in the thermodynamic 
limit, namely the OP is not self-averaging. 

Turning now to G, it is known^iiSS that in the SK 
model the following relation holds for T < T^. 

lim G{T,L) = 1/3. (6) 

L — >oo 
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According to Guerra'^", this relation should hold (for 
T < Tc) in any model which is "stochastically stable" 
with respect to a mean-field perturbation and which has 
a non-self-averaging OP. Under the hypothesis (z) and 
(ii) above, the more general conjecture has also been 
made~° that the above relation holds for T < Tc even if 
the OP is self-averaging. In this case, G would be finite 
but both the numerator and the denominator in Eq.Q 
would vanish, as for example in the Migdal-Kadanoff spin 
glass (see Bokil et al. in Ref. llGj) and the SK spherical 
model""^. It has also been explicitly proven^, under the 
hypothesis {i) and (ii), that one has G{T = Q,L) ~ 1/3 
for any L. Note that models in which G{T = 0, L) ^ 1/3 
in general will not satisfy the conjecture of Ref. jZfl. 



III. NUMERICAL RESULTS 

We simulated the various models with the parallel tem- 
pering technique^, which allows to reach significantly 
lower temperatures than conventional Monte Carlo meth- 
ods. The parameters of the simulation are given in Ta- 
ble ^ Equilibration of the Monte Carlo runs was tested 
by monitoring all the measured observables on a loga- 
rithmic time scale, checking that they had all converged 
within their statistical errors, and by applying the equi- 
libration test discussed in Ref. l26l 



A. Parameter A 

In Figures ^ |21 and |21 we show our numerical results 
for A in the 3D NN, 4D NN, and 3D NNN models, re- 
spectively. The vertical lines in Figures ^ and |21 indi- 
cate the estimated value of Tc. In all cases, the be- 
havior of A resem bles that observed in the SK model 
(see Refs. Il4ll8ll9() . At high temperatures, A decreases 
with L, approximately as , showing that the OP is 
strongly self-averaging in this regime, as expected. Near 



TABLE I: Parameters of the simulations. L is the linear 
system size, Tj^iji and Tmax the smallest and largest tem- 
peratures considered, Nt the number of temperatures in the 
parallel tempering algorithm, A'^s the number of independent 
realizations of the disorder (samples), and MCS the number 
of Monte Carlo steps per spin and per temperature. 



T — Tc, there is a maximum whose position shifts to- 
wards Tc as L increases, an effect of finite-size corrections. 
The shift is modest in the 4D NN model but quite no- 
ticeable in the 3D NN model, where even for the largest 
L the position of the maximum is still significantly larger 
than Tc- In the 3D NNN model, the position of the max- 
imum is also larger than Tc (see discussion in Sec. IIII CI 
on the value of Tc in this model), but the shift is less 
pronounced. The height of the maximum increases with 
L in all models, indicating that A attains a finite value 
in the thermodynamic limit (since it is bounded from 
above), namely that the OP is not self- averaging at Tc, 
as expected. 

At low temperatures (T < Tc), FiguresnEl and|31show 
that A is approximately linear in T, in agreement with 
Eq. (O. Most interestingly, the data for different values 
of L tend to superimpose to each other. In a scenario 
with 6' > 0, the data should tend to zero for large L in 
the whole region below Tc- In 3D, we see no decrease at 
all in the data with increasing L, while a modest decrease 
is observed in 4D. 

To analyze in more detail the size dependence of A at 
low temperatures, in Figure 0] we plot the ratio A/T as 
function of L at different temperatures for the three mod- 
els. The straight lines represent the scaling law Eq.ljSJ as- 
suming 6' ~ 9, and using the estimates of 9 from domain- 
wall calculations, 9 = 0.2 in 3D4-5-7 and 9 = 0.7 in 4D'^. 
No estimates of 9 are available for the 3D NNN model, 
so we use that for the NN model (we expect that is a 
universal exponent equal for both models). Clearly, the 
data in Figure 01 do not agree with the hypothesis 9' = 9 
for the range of sizes considered, and seem to saturate 



4 




0.5 1 1.5 2 

T 



0.2 







1 1 1 1 


I 1 1 1 1 1 1 1 1 


3D NNN - 


- 
- 
- 


L 


i\m 
\ \ \ V\ " 




□ 4 






* 6 






^ 8 






■ 12 




1 1 1 1 


* 16 

1 1 1 1 1 1 1 1 





2 2.5 3 3.5 4 4.5 5 
T 



FIG. 1: Parameter A for the 3D NN model as a function 
of the temperature, for different system sizes L. The vertical 
line represents the estimated value of the critical temperature, 
= 0.95 ± 0.04. 
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FIG. 2: Same as Figure [H but for the 4D NN model, 
vertical line corresponds to Tc = 1.80 ± 0.03. 
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to a constant value instead. We fitted the data with the 
form A(r, V)IT = aL^^ , where 9' should be seen as an 
"effective" exponent which depends on the temperature 
and which effectively takes into account corrections to 
the leading scaling behavior, with 9' — + 0' in the limit 
TL"^' 0. The fits give 9' varying from 0.03 ± 0.02 
(T = 0.7) to 0.00 ± 0.06 (T = 0.2) in the 3D NN model, 
from 0.30 ± 0.05 (T = 1.0) to 0.003 ± 0.006 (T = 0.32) 
in the 4D NN model, and from 0.03 ± 0.04 (T = 2.8) to 
0.08 ± 0.04 (T = 2.0) in the 3D NNN model. 

Therefore, in all cases the data is compatible with 
9' = 0, in agreement with the "many-state" picture, and 



FIG. 3: Same as Figure[T]but for the 3D NNN model. 



is statistically incompatible with with 9' — 9, in disagree- 
ment with the "two-state" picture. As usual, we cannot 
exclude a crossover— to a larger value of 9' for larger L. 
In this case, in the large-volume limit A would be zero 
at all temperatures, except at T = Tc- 

A value of 9' compatible with zero was also obtained 
from the OP distributioni2*2SiSS422i2& and from direct 
measurements of the energy of low-lying excitations cre- 
ated by perturbing the ground state2*i2i. 
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FIG. 4: Log-log plot of A/T versus L at different tempera- 
tures for the 3D NN model (top), the 4D NN model (middle), 
and the 3D NNN model (bottom). The straight lines of slope 
-0.2 (top and bottom) and -0.7 (middle) are the expected 
scaling behavior of Eq.JSJ if 6' = 9. 
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B. Parameters G and B 

Figures|Slin|and[7|show our numerical results for G and 
B for the three models. At high temperatures, G and B 
vanish approximately as in all three models, again 

indicating strong self-averaging. Near T = 0, the data for 
the NN model in both 3D and 4D arc compatible with 
G{T = 0, L) = 1/3 and B{T = 0, L) = 1, as expected^Si^ 
for a continuous coupling distribution. More importantly, 
G{T,L) seems to converge to the value 1/3 for large L 
in the whole low-temperature region, in agreement with 
Eq.®. This is particularly evident in the 4D NN model 
(Figure O, where G{T,L) has already converged to 1/3 
for i = 5 at temperatures below T « Tc/2 « 0.9 (the 
data points above 1/3 are due to incomplete equilibra- 
tioniS.) . It is quite unlikely that the saturation is a finite- 
size effect, therefore our results strongly suggest that in- 
deed G = 1/3 in the whole spin-glass phase. As discussed 
above, if A remains finite as L —^ oo (as indicated by our 
data), this is an expected consequence of stochastic sta- 
bility . If A vanishes, instead, our results would support 
the more general conjecture of Ref. liol 



C. Critical region 

In this subsection we comment on the behavior of G 
and B near the critical temperature, starting from the 
NN model. 

In 3D, the vertical lines in Figure [S] indicates the po- 
sition of the critical temperature, using the value Tc — 
0.95±0.04 quoted in Ref.jS^ which was obtained from the 
parameter B measured in a large-scale simulation. One 
sees that the data for both B and G for different values 
of L come together as T approaches Tc from above, as 
indicative of a phase transition. Below Tc, the data for 
G separate again in a statistically significant way, while 
for B one would need a substantially larger statistics (or 
larger sizes) to see a clear separation, as observed in pre- 
vious studies^^*^. For example, at T = 0.82 the sep- 
aration between the L = 12 and the T = 4 data is 1.4 
standard deviations for B and 2.5 standard deviations for 
G. The small separation below Tc is probably due to the 
vicinity of D = 3 to the lower critical dimensionSiS*2i. 
We also note that the crossing point of G is at temper- 
atures larger than Tc, and close inspection shows that it 
shifts towards Tc from above as L increases. A similar 
shift was observed for the position of the maximum of A 
in Figure n 

In 4D, both B and G display a very clear crossing 
(see Figure EJ, as also observed in previous studiesS2i^. 
From B we estima te Tc = 1.80 ± 0.03 in agreement with 
the results of Refs. I34l42l This value is indicated by the 
vertical lines in Figure El As in 3D, the crossing point of 
G is at temperatures larger than Tc and shifts towards 
Tc as L increases. 

Overall, this confirms that both in 3D and 4D the cor- 
rections to scaling are significantly larger for G and A 



than for B. Since G and A have also much larger statisti- 
cal errors than B, the latter quantity is to be preferred to 
G and A for locating Tc in models with TRS. As already 
mentioned, a much more accurate quantity for this pur- 
pose is the correlation length, which shows a very clear 
crossing in 3D^4, unlike B and G. 

Finally, in the 3D NNN model both G and B show 
a rather clear crossing (see Figure O . This provides a 
clear evidence for the existence of a phase transition in 
3D Ising spin glasses, confirming recent results for the NN 
modcl^''-'^^-^- that obtained a convincing evidence (espe- 
cially Ref 24) after the issue had remained unsolved for a 
long time. The crossing point is at T ~ 3.3 for B and at 
somewhat higher temperatures for G, although also here 
the crossing for G shifts to the left as L increases. From 
the data for B one might be tempted to conclude that the 
critical temperature is Tc ~ 3.3. However, if this was the 
case, the value of B at Tc (which is a universal quantity) 
would be lower in the 3D NNN model than in the 3D 
NN model, violating universality. This suggests that the 
actual value of Tc is significantly lower than 3.3, despite 
the clear crossing of B (which would then be strongly 
affected by scaling corrections), and for this reason we 
have not indicated the position of Tc in Figures 13 and |7| 
A more detailed analysis^ clearly shows that indeed Tc 
is significantly lower than 3.3 in this model. 
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FIG. 5: Parameters G (main figure) and B (inset) for the 3D 
NN model, as a function of the temperature and for various 
system sizes. The vertical lines correspond to Tc = 0.95±0.04, 
the horizontal line in the main figure corresponds to the limit 
relation G = 1/3. 



IV. CONCLUSIONS 

To conclude, we have provided evidence that the order 
parameter is not self-averaging in the low-temperature 
phase of the Edwards- Anderson Ising spin glass in 3D 
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FIG. 6: Same as Figure \^ but for the 4D NN model. The 
vertical lines correspond to Tc — 1.80 ± 0.03. 



and 4D, which implies an exponent 9' ~ 0, in agreement 
with a "many-state" picture of the spin-glass phase, such 
as the repUca-symmetry-breaking picture or the "TNT" 
picture. As usual, due to the limited system sizes that are 
currently reachable in numerical simulations, we cannot 
exclude that for larger sizes one recovers self-averaging, 
nevertheless our result is consistent with other studies 
which used sample-averaged quantitie3242i2S42i and also 
found 6*' ~ 0. Independently of whether there is self- 
averaging or not, we have provided evidence that the 
identity G = 1/3 holds in the thermodynamic limit in 
the whole spin-glass phase, a fact that calls for a theo- 
retical explanation in terms of the geometry and ener- 
getics of the low-lying excitations. We have confirmed 
that G and A can be used to locate the spin-glass tran- 
sition, although in models with time-reversal symmetry 
the usual sample-averaged parameters provide a better 
determination. Finally,we have confirmed the existence 
of a spin-glass phase transition at finite temperature in 
three dimensions. 
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FIG. 7: Same as Figure but for the 3D NNN model. The 
true critical temperature is significantly lower than the cross- 
ing point of B. 
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